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I .   INTRODUCTION 

A.   BACKGROUND 

Extremely  low  frequency  (ELF)  geoelectric  noise  has  a 
highly  impulsive  nature  which  makes  it  very  difficult  to 
remove  by  filtering.  The  main  contributors  to  the 
impulsiveness  of  the  ELF  spectrum  are  lightning  discharges. 
The  lowest  frequencies  of  the  ELF  range  (1  to  60Hz)  have  such 
a  low  attenuation  rate  as  they  propagate,  that  lighting 
discharges  all  over  the  world  must  be  considered  when  studying 
the  noise  spectrum  at  a  single  location.  ELF  signals  have 
been  studied  for  such  applications  as  power  transfer, 
communications,  geophysical  surveys,  and  so  forth.  The  amount 
of  information  on  the  subject  is  immense  and  increasing.  How 
to  account  for  the  impulsive  background  electromagnetic  fields 
inevitably  present  in  any  system  operating  in  this  frequency 
range  remains  a  serious  problem  that  any  system  designer  must 
answer . 

The  background  noise  of  the  lower  ELF  frequency  range  can 

be  thought  of  as  consisting  of  both  local  and  global 

phenomena.   The  global  component  is  comprised  of  signals  that 

are  relatively  coherent  over  extended  distances.   A  detailed 

recording  of  data  in  this  frequency  range 

...is  characterized  by  bursts  of  a  few  cycles  duration 
which  may  be  ten  or  more  times  the  mean  amplitude  between 


bursts.  This  bursty  character  is  apparently  coherent  over 
most,  if  not  all,  of  the  Earth. [Ref .  1] 

The  local  phenomenon  is  comprised  of  those  components  that  are 

not  coherent  over  extended  distances  and  are  assumed  to  be 

products  of  the  local  environment  of  the  detector.   A  set  of 

parallel   detectors   separated   by   some   distance   can   be 

coherently  combined  to  remove  the  global  portion  of  the 

background  noise  and  thus  lower  the  noise  floor. 

B .  OBJECTIVE 

The  Applied  Physics  Lab  (APL)  at  Johns  Hopkins  University 
has  for  years  studied  and  made  measurements  of  the  background 
geoelectric  noise  on  widely  separated  pairs  of  electrodes. 
The  purpose  of  this  thesis  is  to  describe  techniques  of 
adaptive  filtering  that  the  author  applied  to  data  collected 
by  APL,  in  order  to  determine  the  amount  of  noise 
cancellation.  The  techniques  are  implemented  using  software 
developed  during  an  extended  visit  to  the  Applied  Physics  Lab 
by  the  author.  In  a  related  effort,  software  will  also  be 
presented  implementing  appropriate  detection  methodology  for 
the  detection  of  narrowband  signals  embedded  in  the  noise. 

C .  ORGANIZATION 

The  thesis  begins  with  a  brief  discussion  of  the 
constituents  of  this  ELF  noise  to  demonstrate  the  difficulty 
in  modeling  of  the  spectrum  and  to  provide  supporting 
background  information.   The  basic  concepts  of  the  adaptive 


filtering  will  be  presented  as  a  lead  in  to  the  adaptive 
algorithm  to  be  used  to  cancel  the  noise.  These  introductory 
sections  will  be  followed  by  a  discussion  of  stability  and 
convergence  concerns  for  the  adaptive  algorithm  in  this  type 
of  data  environment.  The  detection  discussion  is  then 
presented  as  a  continuation  of  the  processing  required  for  the 
detection  of  narrowband  signals.  The  adaptive  algorithm  is 
then  used  in  the  detector  scheme  to  demonstrate  performance  of 
the  developed  software.  A  final  section  will  summarize  the 
work  performed,  draw  relevant  conclusions  about  the  work  and 
propose  possible  future  studies  that  could  be  derived  from 
this  work. 


II.  THE  ELF  BACKGROUND 

Many  studies  have  been  conducted  in  an  attempt  to  develop 
an  understanding  of  what  actually  constitutes  the  natural  and 
man-made  background  noise  spectrum  in  the  one  to  sixty  hertz 
range.  The  sources  of  the  energy  in  this  frequency  range  can 
be  primarily  attributed  to  atmospheric  disturbances  and  man- 
made  power  transfer  systems.  A  major  problem  with  these 
sources  of  noise  is  that  they  are  highly  irregular  in  both 
magnitude  and  duration  of  the  noise  signal.  Attempts  to 
classify  and  assign  specific  characteristics  to  these  sources 
have  proven  difficult.  A  brief  history  of  the  research 
conducted  in  these  areas  is  provided  in  the  following  two 
subsections . 

A.   ATMOSPHERIC  NOISE 

The  source  of  the  majority  of  noise  in  the  one  to  sixty 
hertz  frequency  range  is  lightning  discharges.  The  lower  ELF 
range  discussed  here  also  contains  what  is  called  the  Schumann 
Range.  The  Schumann  range  takes  its  name  from  the  early 
theoretical  work  of  the  German  scientist  W.  0.  Schumann. 
Schumann  theorized  that  the  cavity  created  by  the  surface  of 
the  Earth  and  the  Ionosphere  has  naturally  occurring  resonant 
frequencies.  He  took  this  research  further  to  derive  a  set  of 
harmonic  resonant  equations  for  calculating  these  frequencies. 


The  work  of  Schumann  has  been  carried  further  with  the 
actual  measurement  of  the  Schumann  resonance  frequencies  in 
field  experiments.  In  1959  and  I960,  H.  L.  von  Konig 
substantiated  the  presence  of  Schumann  resonances  by  observing 
the  waveforms  in  the  output  of  a  narrow  band  amplifier.  The 
experimental  presence  of  Schumann  resonances  created  a  large 
quantity  of  literature  on  the  subject  as  scientists  attempted 
to  understand  the  phenomenon.  Once  a  foundation  of  the 
propagation  properties  of  the  Earth-Ionosphere  cavity  was 
established,  the  search  for  the  source  of  the  excitation 
began. [Ref .  2] 

The  most  prominent  excitation  source  of  the  Earth- 
Ionosphere  cavity  is  the  cloud  to  ground  lightning  stroke. 
The  number  of  similar  strokes  present  in  any  one  flash 
observed  during  a  thunderstorm  is  a  random  variable  that 
normally  ranges  from  two  to  twelve.  Since  the  number  of 
strokes  is  generally  too  high  to  measure  individual  waveforms, 
a  more  applicable  measure  is  the  power  spectral 
density  [Ref .  3]  .  The  frequency  spectrum  of  this  noise  was 
resolved  in  sufficient  detail  to  estimate  the  quality  factor, 
Q,  of  the  earth-ionosphere  cavity  by  M.  Balser  and  C.  A. 
Wagner  in  1960.  The  Q  of  a  cavity  resonator  is  a  measure  of 
the  bandwidth  of  the  resonator  [Ref .  4].  The  source  of  the 
noise,  as  theorized  by  H.  Raemer  in  1961,  is  the  response  of 
the  earth-ionosphere  cavity  to  electrical  discharges  created 
in  thunderstorms  all  over  the  world.   Raemer  attempted  to 


model  the  ionosphere  in  an  attempt  to  reproduce  mathematically 
the  measured  response  of  Balser  and  Wagner.  Although  Raemer ' s 
attempt  fell  short  of  its  goal,  it  began  the  process  of 
refinement  of  a  working  model  of  the  ionosphere  that  continues 
today. [Ref .  2] 

The  early  work  performed  by  these  scientists  brought  out 
the  difficulty  in  attempting  to  model  this  naturally  occurring 
phenomenon  due  to  its  randomness.  Parameters  must  be 
estimated  by  some  method  in  order  to  model  this  random 
process.  The  highly  variable  nature  of  these  parameters  makes 
the  task  of  modeling  the  noise  nearly  impossible.  Examples  of 
these  varying  parameters  are  diurnal  variations  in  the 
ionosphere,  twenty-four  hour  variations  in  the  ionosphere, 
seasonal  variations  in  the  locations  of  thunderstorm  regions, 
and  the  randomness  of  individual  lightning  strikes.  These  are 
only  a  few  examples  of  the  parameters  thai,  must  be  considered 
in  order  to  model  the  atmospheric  noise. 

B.   MAN-MADE  NOISE 

The  presence  of  atmospheric  noise  is  not  the  only  concern 
that  must  be  addressed  when  discussing  background  noise  in  the 
one  to  sixty  hertz  range.  Many  components  in  today's  world 
emit  unshielded  interference  that  falls  in  this  range.  These 
components  are  the  result  of  either  faulty  equipment  or  poor 


design.  Examples  of  this  noise  source  are  power  lines  that 
are  used  to  distribute  alternating  current  through  populated 
areas . 

Large  ferric  objects  can  create  fields  around  them  that 
can  be  detected.  The  strength  of  these  fields  can  be  well 
below  the  previously  mentioned  sources  of  noise  but  can  be 
discernable  in  certain  situations.  Motion  of  these  large 
magnetic  objects  in  the  Earth's  natural  magnetic  field  is  the 
source  of  a  small  field. 

Motion  of  the  detector  can  create  an  increase  in  the 
background  noise.  Any  slight  motion  of  the  detector  can 
create  misleading  information  and  in  most  cases,  an  increase 
in  the  noise  signal.  The  background  ELF  spectrum  can  be 
thought  of  as  the  summation  of  three  orthoganal  vectors. 
Relative  motion  between  these  three  vectors  and  the  detector 
creates  a  false  signal  that  can  raise  the  noise  floor.  This 
type  of  noise  can  be  compensated  for  by  the  use  of  motion 
sensor . [Ref .  5] 

Other  sources  of  man-made  noise  can  be  found  in  the 
detector  equipment  itself.  Improper  shielding  of  cables  and 
electronic  devices  used  in  the  design  of  the  detector  and  its 
processing  components  can  lead  to  elevated  noise  floors. 
Unlike  naturally  occurring  atmospheric  noise,  these  sources 
can  be  traced  and  eliminated.  The  removal  of  this  noise  can 
come  in  two  forms,  elimination  of  the  source  or  isolation  of 


the  source  from  the  detector.  The  fact  that  man-made  noise 
sources  exist  must  not  be  overlooked  when  attempting  to 
process  signals  at  the  output  of  the  detector. 


III.  ADAPTIVE  FILTERS 

A.   GENERAL  TERMINOLOGY 

The  ELF  background  noise  has  been  described  as  a  highly- 
impulsive,  non-stationary  random  process .  Some  of  the  causes 
of  these  variations  have  been  briefly  described  above.  The 
use  of  traditional  filter  designs  may  not  be  the  optimum 
method  for  reduction  of  the  noise  floor  for  reception  of  low 
strength  signals  in  the  one  to  fifty  hertz  range.  Since  their 
introduction,  the  use  of  adaptive  filters  has  shown  promise  in 
combating  exactly  this  type  of  scenario.  To  understand  the 
application  of  adaptive  filters  in  this  case,  a  brief 
discussion  of  relevant  terms  must  be  conducted. 

An  adaptive  filter  is  by  definition  a  filter  that  has  the 
capability  to  adjust,  by  itself,  a  set  of  design  parameters 
that  are  based  on  estimated  statistical  characteristics  of 
the  signal  to  be  filtered.  This  process  of  adjusting  the 
filter  to  the  present  situation  alleviates  the  need  to  have  a 
priori  knowledge  of  the  relevant  signal  characteristics.  The 
previous  chapter  described  how  difficult  it  is  to  derive  a  set 
of  statistical  characteristics  for  the  ELF  noise  environment. 
If  a  priori  knowledge  of  the  relevant  signal  characteristics 
were  known,  an  optimum  filter  could  be  designed.  The  goal  of 
adaptive  filtering  is  to  estimate  these  statistical  parameters 


and  to  refine  the  estimation  through  the  use  of  an  algorithm. 
If  an  appropriate  algorithm  can  be  found,  then  after  a 
sufficient  number  of  iterations,  the  adaptive  filter  should 
converge  to  the  optimum  filter.  The  use  of  adaptive  filters 
requires  the  designer  to  look  for  the  optimum  algorithm  for 
estimating  the  relevant  signal  characteristics  instead  of 
attempting  to  discover  the  actual  signal 
characteristics . [Ref .  6] 

Adaptive  filters  are  divided  into  two  basic  categories, 
open-loop  and  closed-loop.  An  open-loop  configuration  is 
usually  a  two  stage  process.  The  first  stage  is  used  to 
"learn"  the  statistics  of  the  relevant  signal.  The  results  of 
the  first  stage  are  then  used  in  the  second  stage  to  compute 
the  filter  parameters  required  using  a  nonrecursive  algorithm. 
In  contrast,  the  closed-loop  configuration  uses  oniy  one  stage 
to  develop  the  filter  parameters:  the  relevant  signal 
characteristics  are  not  explicitly  estimated  but  are  inferred 
from  a  recursive  algorithm  that  updates  the  filter  parameters 
directly  as  each  new  data  point  is  obtained.  The  adaptive 
filter  gains  additional  knowledge  of  the  signal 
characteristics  during  each  iteration.  The  gain  in  knowledge 
results  in  an  improvement  of  the  filter  parameters  which  are 
adjusted  and  their  performance  is  used  as  an  input  to  the  next 
iteration  of  the  algorithm.  The  adaptive  filter  discussed  in 
this  thesis  falls  in  the  closed  loop  category.  The  advantage 
of  closed-loop  over  open-loop  is  that  the  requirement  for  only 
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one   stage   generally   corresponds   to   a   less   expensive 
configuration  based  on  its  simpler  implementation. [Ref  6] 

The  description  of  the  filters  that  will  be  analyzed  in 
this  study  carries  terminology  that  is  pertinent  to  all 
filters,  not  just  adaptive  filters.  Filters  are  implemented 
in  either  continuous-time  or  discrete-time  which  defines  the 
form  of  the  relationship  between  the  input  and  output  of  the 
filter.  The  filters  that  will  be  studied  are  of  the  discrete- 
time  version  which  means  that  the  filter  may  be  described  by 
a  difference  equation.  The  structure  of  the  filters  to  be 
studied  are  referred  to  as  finite-impulse  response  (FIR) , 
tapped-delay-line,  or  transversal  filters.  This  description 
means  that  the  filter's  output  relies  only  on  the  past  and 
present  values  of  the  input.  The  advantage  of  a  FIR  filter  is 
its  inherent  stability . [Ref  6] 

The  algorithm  that  is  used  to  update  the  filter  parameters 
is  generally  named  after  some  of  its  operating 
characteristics.  The  algorithm  used  here  is  called  the 
Sequential  Regression  Algorithm  or  SER[Ref.  7].  This 
algorithm  uses  an  iterative  approach  to  approximating  the 
weighting  factors  required  for  a  linear  combiner.  Figure  1 
shows  a  basic  diagram  of  a  linear  combiner. 

This  system  uses  simultaneous  samples  of  the  input  signal 
at  time  index  k,  to  produce  an  output  signal,  yk,  by  a  linear 
operation  with  weighting  factors,  wl .  A  second  way  to 
physically  interpret  the  linear  combiner  is  to  look  at  it  in 
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Figure  1  Basic  Linear  Combiner 

the  form  of  a  transversal  filter.  Figure  2  shows  how  this 
system  is  designed:  X,,  corresponds  to  the  last  L  samples  of 
x  at  time  index  k.   The  weighting  factors  are  dual  indexed  on 


Figure  2  Basic  Transversal  Filter 


the  time  delay,  1,  followed  by  the  time  index,  k.  The  weight 
factor,  w2k,  corresponds  to  the  weight  factor  for  the  sample 
taken  two  samples  ago  from  the  present  sample  time,  k. 
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B.   THEORY  AND  DESIGN 
1.   Configuration 

The  filter  to  be  developed  here  is  a  combination  of 
the  previous  two  figures.  The  transversal  filter  is  applied 
to  a  set  of  multiple  sensor  outputs  to  form  an  input  vector 
that  is  N  x  L  data  points  long,  where  L  is  the  number  of 
separate  inputs  or  sources,  and  N  is  the  length  of  the 
transversal  filter.  There  are  two  reasons  for  using  a  filter 
design  of  this  nature.  The  first  is  the  ability  to  use  more 
than  one  source  of  input,  allowing  for  a  better  spectral 
sampling  of  the  background  ELF  energy.  The  second  is  the 
ability  to  use  the  previous  N  samples  thereby  smoothing  the 
output  since  it  is  not  dependent  on  instantaneous  values  which 
may  vary  greatly  from  sample  to  sample.  Both  of  these  ideas 
are  beneficial  to  the  overall  performance  of  the  filtering 
process . 

The  output  of  the  filter,  yy,  can  be  expressed  in 
terms  of  the  input  matrix  and  the  weighting  factors  matrix. 
To  do  this,  some  defining  equation  must  be  expressed  for  the 
input  matrix,  Xk, 

x*  =  [xlk  x2k    ...  x^V  (l) 

and  weight  factors  matrix,  Wk/ 

***[*!*  w2k  ...  fc^r  (2) 

where  each  sensor's  input  vector  is 
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Xlk  =    [jqdt)   jqdr-1)   x1(/:-2)    ...   xx  (*- (tf-1)  )  ]  (3) 

and  weighting    factor   vector    is 

^1*    =     1^1*       Wl{k-1)        Wl(k-2)         •'•        Vl(ic-(N-1))]    ■  (*) 

The  first  subscript  of  Equations  3  and  4,  labeled  1,  2,  3, 
.  ..,  L,  indicates  the  source  of  the  sample.  The  second 
subscript,  k,  indicates  the  time  of  the  sample.  A  total  of  N 
sample  values  from  each  sensor  are  used  for  each  calculation 
of  yk.  A  physical  model  for  the  filter  to  be  studied  is  shown 
in  Figure  3.  The  scalar  output  value,  yk  in  Figure  3,  is 
defined  as: 

yjt=XjwJt=WjXJt.  (5) 

These  basic  equations  will  become  the  starting  point 
for  the  development  of  an  operating  filter  system  for  the 
removal  of  the  ELF  background  noise.  It  is  preferable  that 
the  weighting  factors,  Wk,  cause  the  output,  yk,  to  be  an 
estimate  of  some  desired  signal,  dy .  The  concept  of  the 
filter  is  to  use  the  weighting  factors  to  model  the 
differences  between  signals  originating  from  the  same  source 
but  sensed  at  various  locations.  The  inputs,  xx  in  Figure  3, 
represent  samples  of  a  signal  taken  at  different  locations. 
The  desired  signal,  dk,  is  also  a  sample  of  the  same 
background  source  signal.  The  difference  between  the  desired 
signal  and  the  output  of  the  filter  is  expressed  as  the  error 
in  the  estimate,  €k. 
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ek  =  <*k  -  Yk 


(6) 


In  the  system  designed  here,  the  estimate  error  is  the 
output  sought  for  analysis.  If  the  weight  factors  are  set 
properly,  they  will  theoretically  remove  all  signals  that  are 
coherent  between  the  different  sample  locations.   The  goal  of 
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Figure  3  Multi-Channel  Transversal  Filter 
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such  a  system  is  to  remove  the  correlated  background  noise 
from  a  single  sensor  thus  enhancing  signals  that  are 
particular  to  that  sensor  alone.  If  the  estimated  value  of  yk 
can  be  made  to  represent  the  ELF  background  noise,  then  the 
estimate  error,  €k,  can  be  assumed  to  be  due  only  to  the 
signal  present  in  the  desired  signal.  The  resulting  output 
from  the  noise  canceler  has  a  much  lower  noise  floor  due  to 
the  removal  of  the  correlated  noise. 
2 .   The  Mean  Square  Error 

The  filtering  process  must  have  a  means  by  which  an 
optimal  solution  can  be  defined.  The  mean-square  error  will 
be  used  to  develop  a  defining  equation,  called  the  performance 
surface,  for  the  optimal  solution.  Combining  Equations  3  and 
4,  a  definition  of  estimated  error  in  terms  of  the  input 
vector,  the  desired  signal,  and  the  weigl  ting  factors  is 
obtained.  The  expected  value  of  the  squared  estimated  error, 
£(ek)2,  assuming  that  €k,  dk,  and  Xk  are  statistically 
stationary  is  given  as: 

E[e2k]    =  E[dk2]  +VTE[XkXTk]V-2E[dkXTk]V  (7) 

The  variables  xk  and  dk  are  assumed  to  be  dependent,  and  in 
this  case  both  contain  the  same  ELF  background.  The  equation 
for  the  mean  square  error,  MSE,  associated  with  the  estimate 
is  thus  defined  in  terms  of  the  auto-correlation  function  of 
the  input  vector  and  the  cross-correlation  function  between 
the  input  vector  and  the  desired  signal. 
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The  auto-correlation  and  the  cross-correlation 
functions  will  be  represented  by  R  and  P  respectively.  The 
elements  of  these  two  matrices  are  constant  values  if  the 
input  vector  xk  and  dk  are  stationary  random  variables.  A 
mathematical  expression  for  the  auto-correlation  matrix,  R, 
is : 


R  =  E  [XA.XJ]  =E 


%lk 


■^2k^: 


X-uptzk 


Ik       ^2k 


Xlk^Lk 

X2k%Lk 


%Lk 


(8) 


XLk^lk    XLk^lk     '  ' 

and  the  cross-correlation  matrix,  P,  is: 

P  =  E[dkXk)   =  S[d^1Jc  d,X2k     .  .  .   d^r  (9) 

Using  these  two  definitions  in  Equation  7  results  in  the 
following  definition  of  the  mean  square  error  in  the  kth 
estimate  of  y: 

MSE  =  J£S  =  E[e^]  =E[dt)  +  WrRW  -  2PrW  (10) 

The  optimum  solution  for  the  selection  of  the  weighting 
factors  is  the  condition  where  the  mean  square  error  is  at  a 
minimum. 

3.   Minimizing  The  Error 

The  minimum  mean  square  error  condition  can  be  solved 
for  by  taking  the  gradient  of  the  MSE  with  respect  to  the 
weighting  factors  and  setting  it  equal  to  zero.   The  use  of 
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the  minimum  value  of  the  MSE  as  the  optimal  value  is  based  on 
the  assumption  that  the  weighting  factors  are  only  designed  to 
model  that  portion  of  the  input  vector  that  is  correlated  with 
the  desired  signal.  If  the  error  is  at  the  minimum  value, 
then  the  filter  is  removing  the  maximum  amount  of  correlated 
data  between  the  input  and  the  desired  signals.  The  system 
being  evaluated  is  based  on  the  assumption  that  only  the 
desired  signal  contains  both  noise  and  a  signal  of  interest 
while  the  input  vector  is  composed  of  a  pure  noise  signal 
alone.  Under  this  assumption,  the  only  correlated  data 
between  the  two  input  data  samples  is  the  noise.  The  desire 
is  to  remove  as  much  noise  as  possible  and  thus  to  remove  the 
maximum  amount  of  correlated  data  between  the  input  vector  and 
the  desired  signal.  It  is  therefore  assumed  that  the  filter 
is  working  at  its  optimum  value  when  the  weiahting  factors 
meet  the  criteria  of  minimizing  the  mean  square  error  of  the 
output . 

4 .   The  Optimum  Weighting  Vector 

The  optimum  weight  factor  will  be  represented  by  the 
symbol  W°.  The  two  conditions  that  must  be  met  in  order  for 
the  value  of  J  to  be  minimized  are: 

VwcJfiS|w=wo  =  0  (11) 


and 
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82  J 

ss    >  0.  (12) 


dwfiWj 


These  two  equations  state  that  when  the  weight  vector  is  at  an 
optimum  the  performance  function  is  at  a  local  minimum  since 
the  slope  is  zero  and  the  curvature  is  upward.  The  gradient 
of  the  performance  function  is  evaluated  to  be: 

V  =  Vw  J=VvE[d^]   -2Vw(WrP)  +Vw(WrRW) 
=  -2P  +  2RW 

Setting  the  gradient  to  zero  and  solving  for  the  optimum 
weighting  factors  results  in: 

W°  =  R_1P  (14) 

This  equation  is  sometimes  referred  to  as  the  Weiner-Hopf 
equation  where  the  weight  factor  is  referred  to  as  the  Weiner 
weight  vector.  The  input  autocorrelation  matrix,  R,  can  be 
inverted  if  it  is  "positive  definite,"  i.e.,  if: 

VrRV  >  0  (15) 

where  V  is  any  non-zero  vector.  When  the  matrix  R  is  "positive 
definite"  then  the  optimum  weight  factors  can  be  solved  for 
directly  if  the  complete  statistical  properties  of  the  auto- 
correlation and  cross-correlation  functions  are  known.  The 
statistical  properties  are  not  known  in  this  case  so  further 
development  is  required  to  obtain  an  estimate  of  the  optimum 
weight  factors . 
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The  second  condition  that  ensures  that  the  error 
estimate  is  at  a  minimum  is  that  the  partial  derivative  of  the 
gradient  with  respect  to  the  weighting  factors  is  greater  than 
zero.  To  prove  this  the  we  take  the  partial  derivative  of 
Equation  3  with  respect  to  W  : 

d2j       -  2R  (16) 


dw^Wj 


This  result  again  shows  that  the  input  autocorrelation  matrix 
must  be  "positive  definite"  in  order  to  solve  for  the  optimum 
weighting  factors. 

The  performance  surface  can  now  be  expressed  in  terms 
of  the  minimum  error  possible,  the  weighting  functions, and  the 
autocorrelation  of  the  input  vector.  The  minimum  mean  square 
error  possible  can  be  found  for  the  case  when  the  weighting 
vector  is  at  its  optimum  value.  The  opt imum  weighting  factor 
is  used  to  solve  for  the  minimum  error  by  substituting  W°  into 
the  minimum  error  equation. 

Jm±n   =  E[dl]   +  [W°]rRW0-2PrW° 

=  E[d£]   +  [R-'P]  rRR-1P-2PTR-1P  (1 

=  E[dl]  -   PrR-aP 
=  E[dl]  -  PrW° 

The  minimum  error  can  now  be  used  as  the  base  line  by  which 
all  other  errors  are  found.  The  error  at  any  time  is  the 
minimum  error  with  the  addition  of  uncertainty  in  the 
components  used  to  derive  the  minimum  error.  Thus: 
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J   =  Jrmin+  (W-W°)rR(W-W°)  (18) 

The  expression  can  be  shown  to  equal  that  of  Equation  10  by 
the  following  proof.  Using  the  fact  that  (A-B)T  in  general 
can  be  expressed  as   AT-BT: 

J  =   Jmin+  [W°]  rRW°+WrRW-WrRW°-  [W°]  rRW        (19) 

since  all  terms  are  scalars,  the  transpose  is  equal  to  the 
scalar  value,  thus  the  last  two  terms  are  equal.  Substituting 
for  the  minimum  error  and  the  optimum  weighting  vector: 

J  =  E[d£]   -P7R"1P+  [W0]rRW°+WrRW-2WrRW0 

=  E[d*]   -PTR-1P  +  P7R-1RR-1P+WrRW-2WrRR-1P       2 
=  E[df]   +WTRW-2WTP 
=  E[df]   +WrRW-2P7,W 

This   result   validates   Equation   18   by   showing   that   it 
equivalent  to  Equation  10.    An  important  note  about  this 
expression  of  the  performance  surface  is  that  it  is  quadratic 
with  respect  to  the  weighting  vector. 
5.   Iterative  Calculation  of  Wk 

The  optimum  weight  factors  can  be  arrived  at 
iteratively  by  use  of  a  gradient  search  of  the  performance 
function.  What  this  means  is  that  after  each  estimate  of  the 
optimum  weight  vector,  WK,  the  slope  or  gradient  of  the 
performance  function,  in  the  direction  of  the  optimum 
weighting  vector,  is  used  to  derive  the  next  weighting  vector, 
Wk+1.   This  method  ensures  that  each  estimation  of  the  optimum 
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weight  factors  is  at  least  as  close  or  closer  to  the  optimum 
value  than  the  previous  estimate.   Multiplying  Equation  13  on 
the  left  by  Vi   R1  and  combining  the  result  with  Equation  14, 
we  obtain: 

W°  =  W-1r-xV  (21) 

2 

To  convert  this  into  an  iterative  process,  the  optimum 
solution  is  considered  to  be  the  value  of  the  weighting  vector 
for  the  next  iteration. 

w*+1  =w*-jR_lv*  (22) 

This  equation  is  referred  to  as  the  Newton  Method  algorithm 
for  determining  the  weighting  vector. 

The  algorithm  takes  only  one  iteration  to  arrive 
theoretically  at  the  optimum  solution  for  the  values  of  the 
weighting  vector.  This  algorithm  would  be  the  ideal  method  to 
use  if  exact  solutions  to  the  gradient  term  and  the  inverse 
correlation  of  the  input  vector  were  known  exactly.  This  is 
not  generally  the  case  and  an  estimate  for  the  gradient  must 
be  used.  To  compensate  for  the  lack  of  knowledge  of  the 
gradient  at  iteration  k,  a  constant  u,  will  be  used  instead  of 
the  factor  H  in  Equation  22.  The  requirements  on  the  range  of 
values  that  u  can  take  on  are: 
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0  <  n  <  1  (23) 

The  effect  of  u  on  the  solution  at  each  iteration  is  that  it 
affects  the  rate  of  convergence  of  the  algorithm.  A 
discussion  of  the  attributes  and  further  limitations  on  the 
value  of  u  will  be  covered  later.  A  general  form  of  the 
Newton's  method  of  gradient  search  is  therefore: 

W*+1  »  W*-/!*-1^  (24) 

This  expression  is  ideal  under  the  following  three  conditions: 

1.  y  =  V2 

2.  Exact  knowledge  of  the  gradient  vector,  Vk. 

3.  Exact  knowledge  of  the  (unchanging)  matrix  R1. 

These  conditions,  unfortunately,  are  not  normally- 
attainable.  The  value  of  u  is  selected  between  0  and  Vz  so 
that  the  algorithm  is  overdamped  and  can  accommodate 
fluctuations  in  the  matrix  R.  The  compensation  for  this 
selection  of  u  is  that  the  algorithm  takes  longer  to  obtain 
its  optimum  solution.  A  second  modification  to  the  value  of 
U  is  to  base  its  value  on  the  average  eigenvalue,  A,av,  so  that 
U  is  replaced  by  u^av 

The  second  condition  pertaining  to  the  exact  solution 
to  the  gradient  of  the  error  surface  is,  in  general,  never 
attainable.  The  gradient  must  therefore  be  estimated  in  some 
manner.   The  estimate  of  the  gradient  that  will  be  used  is 
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based  on  the  square  of  the  present  value  of  the  error  signal 
and  not  the  true  ensemble  average  as  it  has  been  up  to  now. 

^  =  i*l  =  2e-^  (25) 

d«k  d»k 

From  the  definition  of  the  error  signal: 

dz,    _    d[dk-XTkVk] 

to*'  d«k  ~      *k 

which  leads  to  the  following  estimate  of  the  gradient: 

Vk  =   -2e^  (27) 

With  these  two  relaxations  from  the  three  ideal  conditions, 
the  algorithm  is  now  in  a  form  known  as  the  Newton/LMS 
algorithm[Ref  7] : 

W*+1  =*k  +  2pAavR-iXk  (28) 

This  is  an  optimum  iterative  algorithm  that  is  only  complete 
if  full  knowledge  of  the  input  correlation  function  is  known. 
This  is  generally  not  the  case  and  thus  a  final  modification 
must  be  done  to  produce  a  working  estimate  to  this  algorithm. 
The  final  modification  is  an  iterative  estimate  of  the  inverse 
of  the  input  correlation  matrix,  R"1,  and  is  called  the 
Sequential  Regression  (SER)  algorithm. 


IV.  THE  SEQUENTIAL  REGRESSION  (SER)  ALGORITHM 

A .   DEVELOPMENT 

The  Sequential  Regression  (SER)  algorithm  uses  the 
LMS/Newton  algorithm  as  its  basis  in  conjunction  with  a  method 
for  iteratively  estimating  the  value  of  the  inverse  input 
auto-correlation  matrix  to  produce  the  weighting  vector  for 
the  filter  coefficients.  The  development  of  the  algorithm 
comes  from  Adaptive  Signal  Processing  by  Widrow  and 
Stearns [Ref  7] .  The  core  of  the  algorithm  is  the  use  of  an 
estimate  for  the  inverse  of  the  input  auto-correlation  matrix 
which  reduces  computational  load  with  minimal  estimation 
error. 

The  input  auto-correlation  matrix  is 

R=£'[XA.XJ]  (29) 

where  the  subscript  k  goes  over  the  entire  length  of  the 
random  process  Xx .  R  will  be  estimated  iteratively  using  only 
a  finite  or  truncated  sample  of  the  random  process  X.  The 
estimate  of  the  input  auto-correlation  matrix  is  given  by: 


Ri.  -  — Z-  X  j X  j . 

k        ic+1  i=o  J 


(30) 


This  estimate  is  unbiased  when  the  input  variable,  X,,  is  a 
stationary  random  process.  When  X}  is  not  stationary,  it  can 
be  shown  that  the  estimate  is  not  very  accurate.  The  estimate 
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of  RK  at  each  value  of  k  is  equally  influenced  by  the 
previous,  k-1,  estimates  of  the  input  auto-correlation  matrix. 
In  order  that  the  most  recent  estimates  of  Rk  have  more 
influence  over  the  new  estimate,  a  "fading"  memory  term,  a,  is 
introduced  which  reduces  the  significance  of  the  past 
estimates.  The  fading  memory  expression  of  the  input  auto- 
correlation matrix  is  denoted  by  Qk  where: 

Qk  =  g  a*-*Xkxl  (3D 

The  value  of  a  is  chosen,  as  a  rule  of  thumb,  such  that  the 
half  life  of  the  exponential  function  is  equal  to  the  length 
of  stationarity  of  X.  [Ref  7]  .  This  rule  of  thumb  leads  to  the 
following  statements  about  the  value  of  a. 

0,  a   2  (-l/lengtho£  stationarity  of  X)  (32) 

0  <  a  <  1  (33) 

The  value  of  the  fading  memory  factor  over  k  iterations  is: 


1  a  (34) 


and  thus  the  estimate  of  the  input  auto-correlation  matrix 
becomes : 

R„  -  -^-Q,   =  -i^SU  V  o^X.xT.  (35) 


k       i-ak+x    k       i-a*+1  h 
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This  estimate  is  exact  for  the  condition  where  Xk  is  constant. 
If  the  input  variable,  X  ,  is  stationary,  a  is  one,  and  if  the 
limit  of  Equation  35  is  taken  as  a  approaches  one,  the  result 
is  in  agreement  with  Equation  30. 

The  derivation  of  the  SER  algorithm  begins  with  the 
estimate  of  the  input  auto-correlation  matrix.  Using  this 
estimate  in  the  equation  for  the  optimum  weighting  factors 
results  in 

M*  -  **■  <36> 

The  definition  for  the  estimate  of  the  cross-correlation 
vector,  Pk,    is,  similarly  to  Equation  35, 

k       l-ak+1  fa 
Using  Equations  35  and  37  in  36  gives: 

w*  =  E«*-J<*i*i.  (38) 


The  SER  algorithm  attempts  to  estimate  the  next  value  of 
the  weighting  vector,  Wk+1,  based  on  the  present  values  of  %: 
and  Pv ..    To  do  this,  Equation  3  8  becomes: 
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k-\ 

1  =  0 

The  last  line  of  Equation  39  is  the  result  of  the  following 
relationship : 

Qk  =  aQk^XkXTk.  (40) 

The  value  of  dv   can  be  replaced  by  recalling  Equation  5  and  6 
which  state  that 

<**  =  **  +  y*  =  e*+XJ"*  (41) 

thus  Equation  3  9  becomes 


(42) 


=  Q^W*  +  e***. 
Multiplying  on  the  left  by  Qk _1,  gives  an  iterative  method  for 
calculating  the  weighting  factor  vector 

**♦!  =WA  +  Q*«A-  (43) 

The  expression  for  Wk+1  in  Equation  43  is  similar  to  the 
LMS/Newton  method  derived  earlier.  With  this  in  mind,  an 
approximation  for  the  LMS/Newton  algorithm  can  be  made  using 
the  above  derivation.  The  SER  algorithm  used  as  an 
approximation  of  the  LMS/Newton  method  is  given  as: 
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2nWl-g**)  Q-ie  (44) 


'k+l 


l-cc 


The  algorithm  requires  a  method  of  iteratively  solving  for  the 
value  of  Qk_1  in  order  to  be  complete. 

The  method  to  iteratively  solve  for  Qk  l  begins  by  pre- 
multiplying  Equation  40  by  Qk  ",  post-multiplying  by  Qk_i_1,  and 
then  by  Xk,  to  obtain: 


Q'kQkQi1-!**  =  QkVQkiQl1-!**  +  Ql'x^lQk^Xk 
which  reduces  to, 

=  Q-^xAa+xlQi^Xj,  ) 


(45) 


(46) 


The  quantity  in  parentheses  is  a  scalar  and  can  be  divided  out 
while  the  quantity  Xv/Q,,.!"1  is  also  multiplied  on  the  right  to 
give 


Qjc-ix k* kQk-i  _  ~  i_  «.7>.-i 

~ — : _  Vk  Ajt*icUjc-i 

*+XkQk-iXk 


(47) 


The  first  lines  of  Equation  46  and  47  can  be  combined  and  re- 
arranged to  form, 


a;1 


a 


0-i  (Q'^k)  (fc-i**)  T 

Ujc-i  ^ — — — 


a+XjUQ^X*) 


(48) 


An  iterative  method  of  calculating  Qk_1  is  now  available  and  is 
only  based  on  its  previous  value,  Qk  V1,  and  the  present  value 
of  the  input,  Xk .   The  value  in  parentheses  is  used  three 
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times,  it  will  be  calculated  separately  and   represented  by 
S(<  • 

The  SER  algorithm  is  now  completely  derived  as  an 
iterative  method  of  approximating  the  LMS/Newton  method  of 
adaptive  filtering.  The  following  summary  shows  the  process 
by  which  each  new  value  of  the  weighting  factor  vector  is 
calculated. 

„    g,   2  (-l/lengthof  etationarity  of  Z) 

S  =  Ql^k 


Y  =  a+xjs 


Q"  =  iK-1_7ss1 


(49) 


2tadv.e(i-«^)    , 

WJc+l    WJt  +  -^OL *  * 

^"max 

The  values  that  must  be  carried  forward  from  one  iteration  to 
the  next  are  Qyi1,  and  W.  .  These  represent  the  previous  value 
for  the  estimate  of  the  inverse  input  auto-correlation  matrix 
and  the  present  estimate  for  the  weighting  factor  vector. 

The  output  of  the  algorithm  is  the  error,  €k,  between  the 
present  value  of  the  desired  signal  and  estimate  of  the 
desired  signal : 
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OUTPUT  =  e^  =  dk-WkXl  (50) 

This  value   is   saved  after  each  iteration   in  a  vector 
representing  the  uncorrelated  portion  of  the  desired  signal 
and  noise  vector  of  the  primary  channel  with  the  noise  only 
vectors  of  the  reference  channels  at  the  time  index  k, 
Xk. [Ref  7] 

B.   INITIALIZATION 

The  SER  algorithm  assumes  that  the  value  of  the  inverse 
input  auto-correlation  matrix  and  the  weight  factors  vector  is 
known  from  the  previous  step.  At  some  point  in  the  past,  an 
estimate  for  these  terms  must  have  been  made  so  that  a 
starting  point  can  be  established.  The  required  accuracy  of 
the  estimates  depends  on  the  length  of  time  that  the  algorithm 
will  be  operating  and  the  earliest  time  that  accurate  data  is 
required. 

The  driving  force  behind  the  accuracy  of  the  initial 
estimate  is  the  speed  with  which  the  algorithm  can  converge  to 
an  optimal  solution.  The  speed  of  convergence  is  dependent  on 
two  elements  of  the  algorithm.  The  first  element  is  the 
initial  value  of  the  inverse  input  auto-correlation  matrix, 
Qk~x.  The  significance  of  this  initial  element  is  that  the 
further  the  initial  estimate  is  from  the  true  value,  the 
further  the  algorithm  has  to  adapt  in  order  to  approximate  the 
optimum  solution.  The  second  element  is  the  step  size,  2uJlave, 
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which  specifies  how  far  the  algorithm  can  progress  towards  the 
optimum  value  in  one  iteration.  The  step  size  will  be 
discussed  in  the  next  section  but  for  now  it  shall  be  assumed 
to  be  a  constant  in  this  discussion  of  the  initial  conditions. 
The  use  of  a  large  constant  multiplied  by  the  identity 
matrix  as  the  initial  value  of  the  Qv  l  has  been  sufficiently 
argued  by  Lee  [Ref .  8] . 

Qo1  =  Qbl  <51> 

As  discussed  by  Lee,  the  use  of  a  large  constant  for  q0  is  the 
proper  selection  for  the  case  where  the  random  process  is 
stationary.  Widrow  and  Stearns  argue  that  this  same  choice 
for  the  initial  value  of  q0,  in  a  non-stationary  condition,  is 
also  appropriate [Ref  7] .  The  point  is  made  that  if  a  priori 
knowledge  of  the  random  process  allows  for  a  better  estimate 
of  the  value  of  Qo"1,  then  that  information  should  be  used. 
The  method  of  initialization  for  this  study  is  to  use  Equation 
51  with  q0  equivalent  to  100. 

The  initial  values  used  in  the  weighting  factor  vector  are 
assumed  to  be  zero.  This  assumption  is  made  since  knowledge 
of  the  process  before  time  zero  is  unknown.  The  use  of  a  32 
tap  delay  filter  also  allows  for  an  estimate  of  the  initial 
full  weight  factor  vector  but  this  is  done  so  at  the  expense 
of  the  first  32  output  data  points.  This  is  done  by  expanding 
the  number  of  elements  in  the  weight  factor  vector  as  each  of 
the  first  32  points  are  read  into  the  algorithm.    The 
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algorithm  is  applied  only  to  the  data  that  is  available; 
therefore,  the  first  time  that  each  element  of  Wr  is  non-zero 
is  at  time  index  31.  This  corresponds  to  the  time  when  the 
first  32  data  points  are  input  into  the  algorithm.  The  result 
of  this  initialization  is  that  the  first  real  data  point  to  be 
considered  an  output  of  the  SER  method  used  here  is  at  k  =  31. 
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V.   STABILITY  AND  CONVERGENCE 

A.   THE  SINGLE  WEIGHT  CASE 

The  SER  algorithm  contains  values  that  can  be  set  to 
affect  the  rate  at  which  the  algorithm  achieves  the  optimum 
solution.  A  simple  case  where  only  one  weight  factor  is 
concerned  will  be  used  to  demonstrate  the  relationships 
involved  in  the  selection  of  the  convergence  rate  factors. 
The  single  weight  iterative  process  using  a  gradient  search 
method  like  the  SER  algorithm  can  be  represented  by: 

*jm  =  ■rjt  +  HC-VJ  (52) 

The  expression  for  the  single  weight  gradient  at  time 
index  k  is  given  by: 


d[X(w-w°)2] 


(53) 

w=wk 


dw 

=  2\(wk-w°) 

The  value  of  the  input  auto-correlation  matrix,  r00  for  the 
single  weight  case,  is  replaced  by  its  equivalent  eigenvalue, 
k=E[x.l]  .  Combining  Equations  52  and  53,  an  iterative  equation 
for  the  weight  factors  is  found: 
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wk+1  =  wk-2\kX(wk-w°)  (54) 

for  the  single  weight  case.  The  terms  in  Equation  54  can  be 
rearranged  to  combine  like  terms. 

wk  =  w°+{l-2\L\)k(w0-w°)  (55) 

This  expression  contains  the  geometric  ratio  of  successive 
terms  which  is  also  used  to  define  the  stability  requirements. 
Stability  is  guaranteed  if  the  absolute  value  of  the 
geometric  ratio  is  less  than  one,  i.e., 

|l-2|iX|<l  (56) 

The  limits  placed  on  the  value  of  the  u  can  be  expressed  in 
terms  of  the  eigenvalue  of  the  input  auto-correlation  matrix 
X   by  rearranging  Equation  56 

^>|i>0  (57) 

The  optimum  solution  can  be  seen  as  the  solution  of  Equation 
55  if  the  limit  of  wk  as  k  goes  to  infinity.  This  is  because 
the  geometric  ratio  approaches  zero  as  k  approaches  infinity. 
The  rate  at  which  the  algorithm  approaches  the  optimal  value 
is  dependent  on  the  value  of  the  geometric  ratio. 

Figure  4  shows  the  effect  of  varying  values  of  the 
geometric  ratio  on  the  rate  of  convergence.  If  the  value  of 
the  geometric  ratio  can  be  maintained  between  zero  and  one, 
then  the  algorithm  will  always  approach  the  optimal  solution 
from  one  side  and  can  be  considered  as  overdamped.   For  a 
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geometric  ratio  of  zero,  the  algorithm  reaches  the  optimum 
value  in  one  iteration  from  one  side  and  can  be  considered  as 
critically  damped.  If  the  magnitude  of  the  geometric  ratio 
falls  between  zero  and  negative  one  then  the  algorithm  is 
underdamped . 
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Figure  4  Effect  of  Geometric  Ratio  (g)  on  the  Rate  of 
Convergence:  overdamped  0<g<l;  critically  damped  g=0; 
underdamped  -l<g<0 


B.   MULTI -WEIGHT  CONDITION 

The  single-weight  condition  discussed  in  the  previous 
section  can  be  expanded  to  develop  an  understanding  of  the 
multiple  weight  factor  case.  The  same  conditions  that  applied 
to  the  single  weight  case  are  true  in  the  multiple  weight 
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case.  The  difference  is  that  each  weighting  factor  has  an 
associated  eigenvalue  and  thus  the  condition  of  stability  is 
now  dependent  on  a  set  of  simultaneous  equations  that  must 
meet  the  stability  condition. 


WJc-W°+ 


(l-2|lA0) 

(1-2JXXJ 
(l-2jiA2) 


(1-211^) 


(WQ-W°) 


(58) 


The  subscript  on  k  indicates  the  index  value  corresponding  to 
a  specific  weight  value.  There  are  a  total  of  L  weighting 
factors  and  thus  there  are  L  individual  eigenvalues,  A.. 

The  value  of  u  in  each  of  the  simultaneous  equations  is 
held  constant.  Since  u  is  the  same  for  all  equations,  the 
eigenvalue  that  has  the  largest  value  will  impose  the  most 
restrictive  conditions  on  the  selection  of  a  value  for  u.  The 
new  range  of  values  that  u  can  take  on  is: 


>^>0 


(59) 

The  condition  stated  in  Equation  59  gives  the  multi-weight 
stability  requirements  that  must  be  met  to  ensure  that  the 
algorithm  will  converge  to  the  optimal  solution. 

The  rate  at  which  the  algorithm  will  converge  was  shown  in 
Figure  4  to  be  dependent  on  the  value  of  the  geometric  ratio. 
Now  that  the  value  of  u  is  a  function  of  the  maximum 
eigenvalue  of  the  input  auto-correlation  matrix,  R,  it  can  be 
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shown  how  this  affects  the  rate  of  convergence.  By  holding 
the  value  of  u  constant  for  each  weighting  factor,  an 
undesirably  slow  convergence  rate  can  occur  for  widely  varying 
values  of  A,.  As  shown  previously,  the  algorithm  converges  in 
one  iteration  when  p  is  equal  to  1/2A  for  the  single  weight 
case.  Using  that  same  criteria,  the  value  of  the  multi-weight 
constant  u  is  l/2Xmax.  The  use  of  this  criterion  ensures  that 
the  algorithm  falls  into  the  overdamped  category  since  u  is 
always  less  than  or  equal  to  1/2^  for  1=[0,  1,  2,  ...,  L] . 
The  weight  factor,  wit  with  the  smallest  value  of  k  will  have 
the  slowest  convergence  rate  since  its  geometric  ratio  will  be 
farthest  from  zero.  This  is  true  because  the  value  of  2uX 
will  be  at  its  minimum  and  will  cause  the  geometric  ratio  to 
approach  one.  When  this  condition  holds,  the  algorithm  may 
take  excessively  long  periods  of  time  to  converge  to  an 
optimal  solution. 

C.   NON-STATIONARITY  EFFECTS 

The  discussion  up  to  this  point  in  dealing  with  stability 
and  convergence  assumed  that  the  process  under  study  was 
stationary.  Under  stationary  conditions,  the  values  of  the 
eigenvalues  of  R  are  assumed  to  be  constant.  A  second  key 
assumption  to  the  development  is  that  the  values  of  the 
eigenvalues  are  known.  The  lack  of  knowledge  of  the  input 
auto-correlation  matrix,  R,  requires  that  some  estimation  be 
made  for  the  value  of  u.   The  value  of  u  must  be  adjusted  to 


ensure  convergence  of  the  algorithm  in  an  environment  where 
the  eigenvalues  of  the  auto-correlation  matrix  are  varying  at 
some  unknown  rate  from  iteration  to  iteration. 

A  method  for  compensating  for  the  variations  in  the  range 
of  values  that  u  can  assume  is  to  use  an  averaged  eigenvalue, 
A.ave,  for  the  selection  of  the  constant  u  for  each  iteration. 
The  value  of  u  is  now  a  constant  for  a  single  iteration  but  is 
changed  from  iteration  to  iteration  based  on  the  average 
eigenvalue.  Stability  of  the  algorithm  is  guaranteed  if  the 
value  of  the  product,  uA.ave,  is  between  zero  and  one.  It  is 
important  to  note  that  the  expression  uAave  is  now  used  only  to 
replace  the  constant  u  used  previously.  The  idea  is  to  scale 
the  average  eigenvalue  for  that  iteration  by  a  constant,  u, 
and  use  it  as  the  factor  for  controlling  convergence  of  the 
algorithm. 

The  geometric  ratio  which  controls  the  convergence  of  the 
algorithm  has  now  been  altered  to  include  a  term  to  compensate 
for  the  non-stationarity  of  the  input  random  process,  X.  The 
new  expression  for  the  geometric  ratio  is 

|1   -  2^XavreR"1|   <  1  (60) 

The  inverse  input  auto-correlation  matrix,  R"1,  is  included  in 
the  development  of  the  SER  algorithm.  The  convergence 
controlling  term  is  thus  uXave,  and  must  be  less  than  the 
maximum  eigenvalue  of  R  to  ensure  stability  and  convergence. 
A  second  condition  that  can  be  imposed  on  the  product  \ikave   is 
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that  its  value  must  always  be  between  zero  and  one.  Using  the 
fact  that  the  spread  between  maximum  and  minimum  eigenvalues 
can  be  rather  large,  the  product  of  ]ikave  must  be  much  less 
than  one . 

The  ELF  background  noise  processes  have  shown  a  tendency 
to  have  a  wide  range  between  eigenvalues  and  to  slowly  change 
over  time.  Because  of  these  two  conditions,  the  value  of  the 
uAavo  must  be  adjusted  over  time  to  account  for  these 
variations.  The  chosen  method  is  to  estimate  the  average 
eigenvalue  by  taking  the  average  of  the  square  of  the  input 
vector,  Xk .  This  value  is  multiplied  by  the  1/100  of  the 
inverse  of  the  maximum  of  the  squared  input  vector. 

ux     =  i  x  a  <61> 

100  [max  (xj]2         L 

The  use  of  this  estimate  for  convergence  control  in  the 
algorithm  ensures  that  the  requirements  for  stability  are  met 
for  each  iteration.  The  1/100  term  ensures  that  the  maximum 
value  of  u  is  not  exceeded.  It  was  derived  empirically  with 
the  use  of  an  actual  set  of  ELF  background  measurements  and 
may  require  adjustment  in  a  different  environment. 

The  choice  of  using  the  present  input  to  the  algorithm  to 
get  an  estimate  of  the  average  eigenvalue  is  used  for 
computational  ease  with  minimal  added  error.  The  use  of  the 
inverse  input  auto-correlation  matrix  should  provide  a  more 
accurate  measure  of  the  eigenvalues;  That,  however,  requires 
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the  inversion  of  the  matrix.  The  prevention  of  matrix 
inversion  is  the  exact  reason  the  SER  algorithm  is  developed, 
so  to  obtain  the  eigenvalues  in  this  way  defeats  the  reason 
for  the  algorithm's  design,  and  is  therefore  counter- 
productive. The  method  has  proven  adequate  and  reasonably 
accurate  in  empirical  work  up  to  this  point,  but  does  leave 
room  for  improvement  and  further  study. 
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VI .   IMPLEMENTATION 

A.   A  PHYSICAL  DESCRIPTION 

The  Applied  Physics  Lab  has  two  identical  sensor  platforms 
submerged  in  shallow  water.  Each  platform  has  electrode  pairs 
for  detecting  the  geoelectric  ELF  background  noise.  The 
system  has  been  in  operation  for  over  two  years  collecting 
samples  of  the  background  noise  environment. 

Ideally,  the  two  platforms  are  parallel  when  the  platforms 
are  placed  on  the  bottom  of  the  bay.  Figure  5  shows  the  ideal 
alignment  of  the  platforms.  The  two  platforms  are  labeled  as 
the  east,  E,  and  west,  W,  platforms.  Each  platform  has  a  set 
of  orthogonal  electrode  pairs,  labeled  X  and  Y.  These 
electrode  pairs  are  orientated  to  receive  the  two  orthogonal 
source  vectors,  Ex  and  Ey ,  in  the  horizontal  plane.  The 
strongest  geoelectric  field  strengths  are  in  the  horizontal 
plane  due  to  the  shift  in  polarization  as  a  wave  crosses  the 
water-air  boundary  from  vertical  to  horizontal. 

If  the  two  platforms  are  perfectly  aligned  and  *~he  noise 
recorded  on  each  are  identical,  the  noise  in  a  specific 
direction,  X  or  Y,  can  be  removed  by  subtraction.  To  cancel 
the  x-direction  noise  simply  subtract  XW  from  XE.  Because  the 
platforms  are  not  perfectly  parallel  and  the  recorded  noises 
are  not  identical  though  highly  correlated,  this  method  is  not 
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Figure  5  Illustration  of  Sensor  Set-Up 


ideal.  The  actual  platforms  can  be  misaligned  up  to  10 
degrees  from  perfectly  parallel.  This  is  due  to  physical 
constraints  such  as  bottom  contours  and  mounting  technique. 
The  result  of  such  a  misalignment  is  that  both  the  XW  and  YW 
sensors  of  the  west  platform  have  components  that  are  coherent 
with  XE  of  the  east  platform.   If  the  electrode  pairs  in 
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Figure  5  are  thought  of  as  vectors,  then  it  is  easily  seen  how 
misalignment  can  cause  coherency  between  XE,  XW,  and  YW. 

A  set  of  real  data  is  traced  through  the  signal  processing 
system  designed  to  accompany  this  sensor  system.  A  calibrated 
source  signal  is  embedded  in  the  signal  XE .  The  path  that  the 
received  signal  follows  can  be  seen  in  Figure  6.   The  output 
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Figure  6  The  ELF  Signal  Processing  Block  Diagram 
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of  each  sensor  is  sampled  at  128  Hz.  The  data  is  bandpass 
filtered  and  resampled  down  to  32  Hz  for  processing.  The  data 
is  then  passed  through  the  SER  algorithm  where  the  coherent 
portions  of  the  west  platform  are  removed  from  XE .  The  output 
of  the  SER  portion  of  the  processor  is  demodulated  and  passed 
through  a  likelihood  ratio  detector  where  a  decision  can  be 
made  as  to  when  the  embedded  signal  is  present. 

B.   INPUTS  TO  THE  SER  ALGORITHM 

The  SER  algorithm  is  used  to  reduce  the  ELF  background 
noise  detected  in  the  primary  submerged  electric  field  sensor 
by  cancellation  with  the  output  of  the  reference  parallel 
electric  field  sensor.  These  parallel  reference  sensors  are 
also  submerged  and  are,  ideally,  assumed  to  be  located  within 
the  same  ELF  background  noise  as  the  primary  channel .  The 
primary  channel  y  represents  sensor  XE  and  the  reference 
channels,  xl  and  x2 ,  represent  XW  and  YW,  respectively.  The 
samples  are  the  actual  output  of  the  sensor  system  used  by 
APL .  The  samples  are  taken  at  the  same  time  without  any  known 
signals  present.  The  sampling  rate  of  the  data  input  for  the 
SER  algorithm  is  32  Hz  and  the  total  length  of  data  is  113760 
points.  These  data  samples  correspond  to  almost  an  hour  of 
background  noise  measurements. 

The  filter  is  a  32  point  transversal  filter  that  is 
applied  to  each  channel  simultaneously.  The  length  of  the 
transversal  filter  means  that  the  previous  31  samples  of  the 
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sensor  are  used  with  the  present  sample  value  for  each 
iteration  of  the  algorithm.  The  goal  of  the  algorithm  is  to 
take  the  output  of  the  two  reference  channels,  xl  and  x2 ,  and 
to  generate  an  estimate  of  the  noise  in  the  primary  channel  y. 
Figure  7  shows  what  the  primary  channel's  data  looks  like  in 
real  time.  The  length  of  stationarity,  used  to  determine  the 
forgetting  factor,  a,  is  assumed  to  be  about  300  data  points 
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Figure  7  The  Primary  Noise  (Y)  and  Calibrated  Source  Signal 

long.   This  value  is  used  based  on  empirical  analysis  of  the 
real  data. 

The  effectiveness  of  the  algorithm  is  analyzed  by  running 
a  set  of  scenarios  with  a  known  signal  inserting  in  the  data 
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samples.  The  scenarios  consist  of  manually  inserting  a  32  Hz 
signal  on  a  5  Hz  carrier  that  is  8192  data  points  long. 
Figure  7  shows  the  components  of  the  primary  channel  signal  in 
a  real  time  representation.  The  primary  channel  is  processed 
by  the  SER  algorithm  while  varying  the  strength  and  location 
of  signals  embedded  in  the  background  noise  y. 

C.   SIGNAL  STRENGTH  RESULTS 

The  algorithm  is  first  run  with  the  location  of  the  signal 
within  the  noise  held  constant.  The  strength  of  the  signal  is 
varied  by  dividing  each  data  point  in  the  signal  vector  by  a 
constant.  A  series  of  four  runs  are  conducted,  each 
succeeding  run  has  half  the  signal  strength  of  the  preceding 
run.  The  initial  run  uses  the  signal  shown  in  Figure  7  with 
the  amplitude  divided  by  500.  The  signal  is  indistinguishable 
prior  to  processing  as  seen  in  Figure  8. 

The  signals  for  each  run  are  centered  at  time  19.2,  32  and 
44.8  minutes  into  the  noise  sample  y.  There  is  no 
significance  to  the  positions  chosen  except  that  the  signal 
data  streams  were  not  allowed  to  overlap.  The  objective  of 
this  first  set  of  runs  is  to  demonstrate  how  the  SER  algorithm 
enhances  the  ability  to  distinguish  small  signals  buried  in 
the  noise.  The  data  is  first  processed  in  the  SER  algorithm 
to  remove  the  coherency  between  the  reference  channels  and  the 
primary  channel.  The  output  is  then  demodulated  leaving  only 
the  in-phase  and  quadrature  components  of  the  non-coherent 
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portion  of  the  primary  channel.  An  identical  reference  set  of 
data  is  demodulated  without  passing  it  through  the  SER 
algorithm  so  that  a  visual  and  quantitative  analysis  of  the 
effects  of  the  SER  algorithm  can  be  made.  The  output  of  the 
demodulator  is  shown  in  Figures  9  and  10.  Note  that  the 
envelope  of  the  signal  is  visible  to  nearly  the  strength  of 
l/2000th,  or  the  third  run  made. 

The  output  of  the  demodulator  is  passed  as  input  into  a 
matched  filter  detector  based  on  the  demodulated  signal.  The 
output  of  this  detector  is  shown  in  Figures  11  and  12  with  the 
SER  and  non-SER  output  for  a  given  signal  strength  shown 


Figure  8  Input  Primary  Channel  Noise  And  Signal 
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together.  Numerically,  the  peak  value  out  of  the  detector  is 
three  to  five  times  greater  when  the  SER  algorithm  is  used  as 
compared  to  when  it  is  not  used. 

A  final  numerical  analysis  is  performed  by  using  a  primary 
channel  without  a  signal  installed  to  get  a  measure  of  the 
noise  present  in  the  output.  The  standard  deviation  of  this 
output  of  the  detector  when  this  noise  only  vector  is  used  as 
input  is  used  as  an  approximation  to  the  noise  power. 
Dividing  the  outputs  of  the  detector  in  Figures  11  and  12  by 
this  estimate  of  the  noise  power,  an  estimate  of  the  signal- 
to-noise-ratio  (SNR)  can  be  made.  The  values  of  the  average 
peak  SNR  for  the  three  locations  in  each  run  are  plotted 
versus  signal  strength  in  Figure  13 .  The  SNR  of  the  SER 
applied  data  runs  about  6  dB  higher  than  the  non-SER  applied 
data . 
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Figure  9  Demodulator  Output  of  In-Phase  Component; 
(a)  Signal/500  (b)  Signal/1000  (c)  Signal/2000 
(d)  Signal/4000 
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Figure  10  Demodulator  Output  of  Quadrature  Components; 
(a)  SIGNAL/500  (b)  SIGNAL/1000  (c)  SIGNAL/2000 
(d)  SIGNAL/4000 
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Figure  11  In-Phase  Detector  Outputs:  (a)  Signal/500 
(b)  Signal/1000  (c)  Signal/2000  (d)  Signal/4000 
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Figure  12  Quadrature  Detector  Outputs:  (a)  Signal/500 
(b)  Signal/1000  (c)  Signal/2000  (d)  Signal/4000 
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D.   DETECTOR  THEORY 

The  process  of  adaptive  noise  cancellation  is  performed 
for  the  purpose  of  increasing  the  signal-to-noise  ratio.  A 
section  on  the  narrowband  detection  schemes  that  have  been 
considered  in  conjunction  with  the  noise  canceler  are 
presented  since  they  are  integral  to  the  analysis  here.  Two 
likelihood  ratio  detectors  are  considered  for  the  detection  of 
na  r rowband  signals. 
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The  first  detector  is  designed  with  the  assumption  that 
not  only  is  the  signal  known,  but  its  peak  within  a  given 
interval  is  constant  relative  to  the  start  of  that  interval. 
This  last  assumption  is  to  assure  that  the  demodulated  in- 
phase  and  quadrature  signals  that  are  used  in  the  likelihood 
ratio  detector  match  those  in  the  noise  canceled  and 
demodulated  output  of  the  sensor  system  described  previously. 
In  the  second  detector,  the  latter  assumption  was  dropped  and 
the  design  followed  that  of  an  optimal  envelope  detector.  In 
both  cases  it  was  assumed  that  the  residual  noise  field  was 
gaussian.  This  assumption,  while  not  strictly  true,  seems  to 
yield  quite  reasonable  results . [Ref .  9] 

The  signal  is  generated  by  a  calibrated  source  with  known 
qualities.  The  modelling  equation  for  the  narrowband  signal 
is  given  by: 

s(t)  =  A  (t)  cos  [2nfct  +  8  (t)  ]  (62) 

For  the  first  detector,  the  phase  term,  0,  is  assumed  constant 
within  any  non-overlapping  interval  0  to  T.  The  interval  of 
0  to  T  is  equivalent  to  the  length  of  the  signal.  The  purpose 
of  this  assumption  is  so  that  a  Matched  Filter  can  be  created 
for  the  in-phase  and  quadrature  components  by  using  the 
demodulated  in-phase  and  quadrature  component  of  the  signal. 
Figure  14  shows  the  demodulated  signals  used  for  the  detector. 
The  first  detector  is  a  multi-channel  Matched  Filter  shown 
in  Figure  15.   The  Matched  Filters,  hi   and  hq,    are  determined 
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Figure  14  Demodulated  Signal  for  Match  Filtering 

based  on  the  requirement  that  the  output  signal-to-noise  ratio 
is  maximized  at  time  T.   The  corresponding  equations  are 


/ 


Ru(t-x)   Rlg(t-x) 
Rgi(t-x)   Rqg(t-x) 


hi  (x) 
hq  (x) 


dx 


si  (T-x) 
Sq(T-x) 


(63) 


where  the  quantities  R  denote  the  various  auto  and  cross 
correlations.  This  detector  is  the  one  used  for  all 
processing  performed  in  this  thesis.  The  input  to  the 
detector  is  the  noise  cancelled  output  of  the  SER  portion  of 
the  processor  after  demodulation . [Re f .  10] 

The  second  detector  is  an  optimal  nonlinear  envelope 
detector.   The  description  of  its  principles  and  design  are 
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located  in  the  Appendix.  The  detector  was  not  used  by  the 
author  in  this  work  relating  to  the  analysis  of  the  SER 
algorithm,  but  it  is  presented  as  complementary  work  relating 
to  the  overall  project. 

E.   SIGNAL  LOCATION 

The  signal  strength  analysis  showed  that  the  signal 
divided  by  2000  is  easily  discernable  from  the  background  at 
the  output  of  the  detector.  The  effect  of  varying  the 
location  of  the  signal  will  allow  analysis  to  get  a  better 
feel  for  the  expected  SNR  based  on  a  larger  sample  size  for 
averaging  than  previously  used.  A  second  benefit  of  this 
analysis  is  to  see  the  effect  of  the  initialization  process  on 
the  SNR.  The  signal  will  be  shifted  by  a  quarter  of  its 
length  for  each  run,  or  2048  points.  The  SNR  is  plotted 
versus  the  location  of  the  center  of  the  signal  in  Figure  16. 

The  data  is  run  with  and  without  the  SER  algorithm  to  get 
a  good  feel  for  the  enhancement  achieved  with  the  SER 
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algorithm.   Three  separate  locations  of  the  signal  will  be 

used  for  each  run  to  reduce  calculation  time.   The  use  of 

multiple  signals  causes  no  effect  on  the  SNR  calculation  as 
long  as  the  signals  are  not  allowed  to  overlap. 


VII.   RESULTS  AND  CONCLUSION 

A.   ANALYSIS  OF  RESULTS 

The  use  of  the  SER  algorithm  enhances  the  ability  to 
detect  narrowband  low  power  signals  buried  in  the  background 
noise.  Figure  17  shows  the  spectrum  of  the  noise  sample  from 
sensor  XE  before  and  after  the  use  of  the  SER  algorithm.  The 
frequency  range  is  from  one  to  ten  hertz,  which  is  in  the 
range  of  concern  for  the  JHU  project.  The  spike  seen  at  four 
hertz  (note  the  x-axis  is  a  logarithmic  scale)  is  caused  by 
the  processing  equipment.  This  is  a  man-made  noise  source  and 
will  be  filtered  out  in  the  final  project  design.  The  figure 
illustrates  the  10-15  dB  of  noise  cancellation  possible  by 
using  this  method  of  adaptive  noise  cancellation. 

The  location  of  the  signal  in  comparison  to  the  start  of 
the  adaptive  process  has  been  shown  to  have  very  little  effect 
on  the  signal  enhancement.  The  system  implementation  will 
most  likely  cause  the  algorithm  to  only  require  initialization 
once  a  week  and  operate  constantly  otherwise.  The 
initialization  process  takes  about  100  iterations  or  1.67 
minutes  before  settling  out  but  is  heavily  data  dependent. 
Even  while  settling,  the  performance  of  the  algorithm  is 
adequate  to  receive  most  signals.  For  the  purpose  of  this 
project,  initialization  should  be  of  no  concern. 


59 


■50  - 


-60 


Power  Spectral   Density 


Raw  Data 


After  SER 


V^^Vlw^Aii^fV^ 


'■*v, 


1.0 


3.2 
Frequency   (Hz) 
-axis  is  log  scale 


10.0 


Figure  17  Illustration  of  SER  Algorithm  Cancellation 
Capability,  Power  Spectral  Density  of  1-10  Hz 

B.   CONCLUSIONS  AND  RECOMMENDATIONS 

It  was  hoped  at  the  beginning  of  the  research  that 
information  could  be  gained  on  the  characteristics  of  the 
noise  environment.  The  stationarity  of  the  background  noise 
environment  was  found  to  be  at  least  nine  seconds.  This  is 
believed  to  be  only  a  characteristic  of  the  specific  data  set 
that  is  used.  The  length  of  stationarity  affects  the  value  of 
the  forgetting  factor,  a,  and  is  implied  by  the  use  of  the  SER 
algorithm.  A  better  feel  for  a  good  average  value  to  use  for 
the  length  of  stationarity  will  result  from  the  use  of  more 
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data  sets.   This  analysis  needs  to  be  conducted  and  could  be 
a  source  of  further  work  in  this  area. 

A  second  area  of  work  that  could  develop  from  this 
research  is  a  study  of  the  cancellation  possible  from  a 
vertical  above-surface  sensor  and  a  submerged  horizontal 
sensor.  The  purpose  of  this  would  be  to  differentiate  between 
surface  and  submerged  activity.  This  may  prove  to  be 
impossible  but  it  would  be  interesting  to  see  the  amount  of 
correlation  possible  between  the  two  locations.  This  type  of 
arrangement  could  also  lead  to  a  wider  distance  between 
sensors  which  could  also  increase  the  knowledge  of  the 
background  environment . 

The  purpose  of  the  research  presented  in  this  thesis  is  to 
develop  a  method  to  adaptively  cancel  the  ELF  noise  present  in 
a  submerged  sensor  array.  The  SER  algorithm  appears  to  be  an 
adequate  method  to  accomplish  this  goal  of  noise  cancellation. 
The  cancellation  ratio  in  the  frequency  range  of  interest  is, 
on  average,  12.5  dB.  This  cancellation  ratio  along  with  the 
use  of  optimal  detectors  will  assist  the  system  in  detecting 
narrowband  signals  buried  deep  in  the  background  noise. 
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APPENDIX 

The  second  detector  used  in  conjunction  with  this  project 
of  measuring  the  geoelectric  background  noise  is  an  optimal 
nonlinear  envelope  detector.  The  symbolism  used  is  slightly 
different  but  efforts  are  made  to  define  all  relevant  terms. 
The  in-phase  and  quadrature  data  (after  noise  cancellation) 
are  denoted  by  r£  and  r£,  respectively.  Similarly  the  in-phase 
and  quadrature  noises  are  represented  by  n£  and  n£.  Assuming 
a  joint  Gaussian  distribution  we  have  the  following  two 
hypothesis : 


q 


H0    :     (  jr1  ,    x*    )     =    (    n1    ,    ri 

Hx    :  (  r1    ,    x*    )  =  (  n1    ,    nq    )  +  (  s1    ,  sq  ) 
where  H0    denotes  the  null  hypothesis  (noise  only)  and  Hl 
denotes  the  alternative  (signal  and  noise  present).    The 
vector  notation  means 

r  =  [ri,r*]    =   [{r^z?) (ri ,  x§)  ] 

where  N  denotes  the  length  of  the  interval  for  which  the 
decision  is  being  made.  It  is  then  clear  that  the  probability 
density  function  for  the  received  data  is  as  follows: 


p  {t\hx)   «  exp 


1  (r-8)TR-1(r-8) 

2  ~nn 


where  Er,r,  denotes  the  covariance  matrix  of  the  noise  vector. 
In  this  particular  case,  it  is  assumed  that  the  in-phase  and 
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the  quadrature  components  have  been  decorrelated  and  whitened 
with  the  same  variance  o2  so 


P  (rl      -_MNexpj 


N 


£  (ri-s^i+T  {zg-ag) 


2na2J  [   2o2 

Introducing  polar  variables  (  r£  ,  #.   )  and  {   it  ,    4%     )  it  is 
found  that 


p{z\Hx)    = 


a 


2710 


2N 


exp 


2o- 


£r,2+£s,2-2£r^cos(<tg 


ir=l      ic=l 


Jt=l 


where  rk2  =  rj.2  +  if-   and  tan  (<{>,)  =   r£  /  r£  .   Integrating 
the  above  over  the  angular  variables  will  give 


P(z\Hx) 


_B 


-r  I  ^A^JC 


2KO 


2N 


expl-^Et^-^Il^hr 


where  J0  is  the  modified  bessel  function  of  order  zero.  The 
vector  r  on  the  left  hand  side  now  only  refers  to  the  various 
envelopes,  r  =  [  rlt  .  .  .,  rN]  .  The  same  is  true  for  the  signal 
envelope.   The  likelihood  ratio  detector  is  found  to  be: 


.  pfrltfj 
X= — .  |  .  =exp 
P(z\h0) 


N         N 

5=£'J  S 


r       S^I 
°l     2~ 


The  log-likelihood  ratio  is  then 


2o2m      &*       I  \  o2  ;j 
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The  first  term  is  a  constant,  dependent  only  on  the  known 
signal,  and  can  be  dropped.  Finally  the  test  statistic  or 
detector  output  is [Ref  10]: 


logA  =  £  log  |j( 


jt=i 


Skrk 
0|     2 

o2 


The  implementation  and  analysis  of  this  detector  scheme  is 
left  as  a  possible  follow-on  thesis  topic. 
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